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Abstract. A database on shape evolution of direct iV-body 
models formed out of cold, dissipationless collapse is generated 
using GRAPE and HARP special purpose computers. Such 
models are important to understand the formation of elliptical 
galaxies. Three dynamically distinct phases of shape evolution 
were found, first a fast dynamical collapse which gives rise to 
the radial orbit instability (ROI) and generates at its end the 
maximal triaxiality of the system. Subsequently, two phases of 
violent and two-body shape relaxation occur, which drive the 
system first towards axisymmetry, finally to spherical symme- 
try (the final state, however, is still much more concentrated 
than the initial model). In a sequence of models the influence 
of numerical and physical parameters, like particle number, 
softening, initial virial ratio, timestep choice, different iV-body 
codes, are examined. We find that an improper combination 
of softening and particle number can produce erroneous re- 
sults. Selected models were evolved on the secular timescale 
until they became spherically symmetric again. The secular 
shape relaxation time scale is shown to agree very well with 
the two-body relaxation time, if softening is properly taken 
into account for the latter. Finally, we argue, that the interme- 
diate phase of violent shape relaxation after collapse is induced 
by strong core oscillations in the centre, which cause potential 
fluctuations, dampening out the triaxiality. 
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1. Introduction 

Formation scenarios of elliptical galaxies are either based on 
merging of galaxies or on violently collapsing stellar systems (or 
on combinations of both). Especially the field ellipticals most 
likely pass through a phase of cold dissipationless collapse. It 
only occurs if the primordial gas cloud forming a galaxy has a 
very efficient phase of star formation at the beginning by which 
practically all the gas is transformed into stars. Thereafter the 
system is not supported against further collapse until a dynam- 
ical equilibrium is reached. The structure of the final state is 
restricted by the time-independent tensor-virial theorem (see 
Binney & Tremaine 1987), which includes the possibility of 
anisotropy-supported triaxial configurations. 
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Merritt & Aguilar (1985) and Min & Choi (1989) have 
shown by A-body simulations that the initially spherical sys- 
tem does not remain spherical during collapse. If approxi- 
mately 2T/\U\ < 0.1, where T and U are the kinetic and poten- 
tial energy of the initial state, the system is unstable against 
the formation of a nearly prolate shape. Aguilar & Merritt 
(1990) performed an exhaustive study of the final shape of the 
prolate configuration as a function of initial parameters as ki- 
netic energy and initial triaxiality. Cannizzo & Hollister (1992) 
varied in a similar parameter study the initial density profile 
and the softening parameter (three values). 

This instability, commonly known as "radial orbit instabil- 
ity" (henceforth ROI) was long believed to be connected with 
some initial amount of anisotropy in the velocity dispersions of 
the system, which either stems from the initial conditions or 
develops during collapse. However, Udry (1993) showed that 
even a very small initial anisotropy suffices to develop ROI, 
and more recently Hozumi, Fujiwara & Kan-ya (1996) pointed 
out that the onset of ROI happens during the collapse of a uni- 
form sphere while it remains completely isotropic. So, although 
it is usually stated that the onset of ROI can be understood as 
a type of Jeans instability perpendicular to the radial move- 
ment of collapse, a complete understanding of its mechanism 
is still lacking. It is interesting to note, however, that some 
remarks made by Bettwieser & Spurzem (1986) and Hensler 
et al. (1995) can be applied here. They showed that homoge- 
neous systems in self-similar collapse remain exactly isotropic, 
whereas in other cases anisotropy develops proportional to the 
gradient of u/r, where u is the radial bulk mass streaming 
velocity. 

ROI is related to the concept of global instability of stellar 
systems, although details of this relation are still unclear. The 
global instability of a spherical system of purely radially mov- 
ing stars had been proven by Antonov (1961) and Polyachenko 
& Shukhman (1981) analytically. Most other analytic or semi- 
analytic stability criteria rely on an energy principle, but can 
be applied only to systems where the distribution function is a 
function of energy alone, i.e. isotropic systems (Antonov 1961, 
1962, Lynden-Bell 1964, Ipser & Thorne 1968, Kulsrud & Mark 
1970). For most anisotropic and more realistic cases, however, 
no rigorous proof of sufficient or necessary criteria for stability 
is known, and one has to rely on numerical simulations, as in 
the papers cited in the previous paragraphs. 

Lynden-Bell (1967) derived by using principles of statis- 
tical mechanics a most probable distribution function as the 



2 



end-state of cold, dissipationless collapse, the Lynden-Bell dis- 
tribution. The assumed nature of the process led to the notion 
of violent relaxation. However, already very early iV-body sim- 
ulations (Cuperman, Goldstein & Lecar 1969) were not able to 
show that such a distribution occurs. Burkert (1990) tried to 
relax the problem by assuming that the relaxation is incom- 
plete, i.e. not all elements of phase space are accessible with 
equal probability, but the concept of violent relaxation thereby 
became not so attractive anymore. There was a remarkable se- 
ries of papers (Wiechen, Ziegler & Schindler 1988, Ziegler & 
Wiechen 1989, 1990), who derived the final state as a function 
of the initial state by an energy principle. However, it has not 
yet been practically applied, as far as the authors of this paper 
are aware of. 

Hence direct iV-body simulations still are the appropriate 
tool to examine cold dissipationless collapse. In recent years 
significant advances in hardware and software for direct N- 
body simulations have been made. On one hand TREE codes 
(Barnes & Hut 1986) have been ported to parallel comput- 
ers and coupled with particle based hydrodynamic simulations 
(TREE-SPH, Dave et al. 1997). On the other hand completely 
different approaches like grid-based schemes (Udry 1993) or 
the self-consistent-field method (SCF) by Hernquist & Ostriker 
(1992) have been succesfully used for stellardynamical appli- 
cations (e.g. in the case of cold collapse Hozumi & Hernquist 
1995). Although e.g. the latter appear to be ideally suited for 
collisionless systems, because they do not use particle-particle 
forces to integrate the orbits, they suffer from their inflexi- 
bility in adaptation to different applications and resolutions 
(basically for every new application an appropriate set of ba- 
sis functions to evaluate the potential has to be searched, and 
the attempt to increase small-scale resolution, e.g. in angular 
direction, increases the CPU time nearly as inhibitively as for 
direct iV-body simulations, since CPU time goes roughly with 
n 2 , where n is the number of terms in the series evaluation for 
tangential resolution, much like a force calculation scales with 
iV 2 where N is the total particle number). 

Thus, even for collisionless systems and applications with 
correspondingly very large particle numbers, direct or at least 
TREE-like iV-body simulations will remain very useful and 
important in the foreseeable future because of their inherent 
Lagrangian nature, yielding high resolution at high density re- 
gions automatically. Moreover, we think it is not fully clear 
what is the physically correct case for a collapsing protogalaxy. 
On one hand, the collapse might occur in a system consist- 
ing of protogalactic clouds of masses as large as say 1O 6 M0 
(cf. e.g. Theis & Hensler 1993, 1995); in such case a collision- 
less approach like in the case of SCF codes does not represent 
the physical situation. On the other hand, even considering a 
collapsing system made out of billions of stars only, some re- 
laxation effects (not necessarily just two-body relaxation) can 
occur on timescales much shorter than the standard two-body 
relaxation time. It is not clear and deserves further comparison 
between the different approaches whether a mean field based 
or a direct iV-body method (using large N and very small 
softening) represent a more accurate representation of the real 
physical system under study. 

For the direct iV-body simulations special purpose comput- 
ers have been very successfully built (Sugimoto et al. 1990, 
Makino et al. 1997) and applied (Makino 1996, Makino & 
Ebisuzaki 1996, Makino 1997). There are different versions of 
such machines e.g. for high precision (collisional stellar dynam- 



ics, HARP, GRAPE-4) and low precision (collisionless systems, 
GRAPE-3, GRAPE-5). We have used the GRAPE 3Af and 
HARP-2 special purpose computers in Kiel and Heidelberg to 
produce a data base on the influence of the softening length and 
the particle number on the evolution of the shape (axis ratios) 
of stellar systems during cold collapse and ROI. Due to the use 
of GRAPE 3Af (4.8 Gflop peak performance) our data could 
cover simultaneously small softening and large particle num- 
ber TV. The use of HARP-2 enabled us to provide data without 
any softening for comparison. Small softening should increase 
two-body relaxation, whereas large N decreases it (compared 
to the dynamical timescale). In a real galactic system e should 
be vanishing for all practical purposes, since the radii of the 
stars are extremely small relative to the mean interparticle dis- 
tance. On the other hand, the particle number is so large as 
well (e.g. 10 11 ) that the standard two-body relaxation time is 
large compared to the age of the universe. Since a direct model 
with N = 10 11 and a vanishing softening length e is impossi- 
ble now and in the near future by a wide margin, generally a 
system with very different N and e is used, very often just in 
the vague hope that e is large enough to render all two-body 
relaxation effects unimportant during the simulation and to 
use N as large as possible with the present computer genera- 
tion (compare Theis 1998). As for the evolution of the shape 
of triaxial galaxies, emerging after cold collapse of an initially 
extended gas-free stellar system, we will show that relaxation 
effects still play an unwanted role in the simulations. As a re- 
sult we conclude that great care should be taken in interpreting 
results of such iV-body simulations as models of real galaxies. 



2. Numerical Results 

2.1. The reference model 

In our reference model (model Al) we started with N = 32768 
particles of total mass M = 1. The particles were distributed 
according to a Plummer density profile 

p(r)=p .(l + (r/r ) 2 r 5/2 . (1) 

The units are normalized to a gravitational constant G = 1, 
a unit mass and a total potential energy U = —1/2. For a 
virialized system the latter choice guarantees a normalization 
to the standard N-body units. However, here we start with 
a cold system (i.e. the total energy E w U) resulting in a 

crossing time t CI = (GM 5/2 )/ y /2\E\' i w 1.03, for a virial ra- 
tio r) vil = 2T/\U\ = 0.04. (T is the total kinetic energy.) The 
mass distribution corresponds for a Plummer model to an ini- 
tial scale radius ro = 37r/16 « 0.589 giving a free-fall time at 
the scale radius of iff (r ) = >/37r/(32Gp(ro)) ~ 0.84. p(r) is 
the average mass density inside a sphere of radius r. The ini- 
tial velocities were chosen according to an isotropic Gaussian 
distribution with a velocity dispersion that results in an initial 
virial coefficient »?vir = 0.04. Thus, the collapsing system is pre- 
pared to form a triaxial system via a ROI (Aguilar & Merritt 
1990). 

The gravitational softening length was set to e = 0.01 
which is a factor of two smaller than the average inter- 
particle distance r n — ro/VN at the center and a factor of 
1000 larger than the average 90° deflection impact parameter 
Po = ro/(2N) « 9- 10~ 6 . In opposition to select e large enough 



Fig. 1. Projection of the particles on the (arbitrary) x—y— plane at different times (from upper left to lower right: t = 0, 1.5, 2.5, 5.0, 20.0, 200. 
The simulation was performed with N = 32768 particles and a softening length e = 0.01 (model Al). 
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to have little undesired two-body relaxation we here want to 
examine the different shape evolution of systems with small 
and large softening after and during the so-called "cold col- 
lapse" . It is our aim by this method to study the influence of 
varying degrees of two-body relaxation on the results. The time 
integration in the simulations was performed with a leap-frog 
scheme using a fixed timestep of At = 0.0015 which is less than 
0.3 per cent of the central free-fall time. This small timestep is 
the price for an acceptable energy error in the numerical simu- 
lation of the collapse of an initially cold configuration when a 
small gravitational softening is applied. The energy conserva- 
tion is in all simulations better than 0.2 per cent, mostly better 
than 0.05 per cent. (Typically the energy error is at maximum 
in the collapse phase, whereas lateron the energy fluctuations 
are an order of magnitude smaller than the largest deviation.) 
The potential was evaluated by a GRAPE 3Af special purpose 
computer which applies the standard Plummer softening, i.e. 
the acceleration of particle i by particle j is calculated by 



Gm, 



(r? +£ 2)3/2 



(2) 



The temporal evolution of the system starts with a strong col- 
lapse. After a typical collapse time tg(r ) of the system the cen- 
tral region begins to reexpand whereas the outer layers are still 
falling to the center because of their larger free-fall time. Three 
collapse times later the inner 75 per cent of the system reach 
an equilibrium state, i.e. the Lagrangian radii remain almost 
constant. The half-mass radius has decreased by a factor of 1.6 
compared to the initial value and the inner Lagrangian radii are 
typically smaller by a factor 3-5. The relatively small decrease 
of the half-mass radius is a result of the change in the shape of 
the system: During the first collapse phase the system remains 
spherical, which is consistent with the results of Hozumi, Fu- 
jiwara & Kan-ya (1996). The shape however changes quickly 
afterwards, when the central region reexpands (Fig. 1). After 
t = 1.79tff (ro) = 1.5 a central bar-like structure can be iden- 
tified which becomes more and more dominant until t = 2.5. 
The shape is now almost fixed in the sense that it does not 
vary significantly on a crossing timescale. However, due to re- 
laxation processes which are discussed below there is a secular 
evolution towards a more spherical mass distribution (Fig. 1). 
The dynamics of the central part can be followed by the evo- 
lution of the core parameters and the density centre. In Fig. 
2a the density weighted core radius r core defined in the same 
way as the density radius in Casertano & Hut (1985) decreases 
strongly during the collapse phase and reaches an equilibrium 
value of 0.04 at t = 2.0. The corresponding core mass drops 
from 0.17 for the initial Plummer model to a final core mass 
M coro = 0.05. Once the equilibrium state is reached, both val- 
ues vary only by 10% during the whole simulated evolution 
(until t = 200) which demonstrates the remarkable constancy 
of the core properties. The motion of the density centre will be 
described later in Sect. 3. 

The flattening of the mass distribution is calculated by 
transformation to the principal axes of the tensor of inertia 
relative to the density centre of the system using different num- 
bers of particles which are sorted by their binding energy. For 
50% of all bound particles one finds that the minor to ma- 
jor axis ratio c/a begins to decrease at t = 0.8 and reaches 
a minimum of 0.36 at t = 3.0 (Fig. 2b). Afterwards the axis 
ratio increases linearily to 0.52 at t = 20. The intermediate-to- 
major axis ratio b/a shows a similar behaviour as c/a: With 



Fig. 2a. Temporal evolution of the core radius (solid) and the core 
mass (dashed) for model Al (N=32768, e = 0.01). 



a short delay of At = 0.2 b/a starts to decrease at t = 1.0 
and reaches a minimum of 0.42 at t = 3.0. Lateron, the ratio 
c/a is on average 0.1 less than b/a which indicates a slightly 
triaxial (but near to prolate) system characterized by a triax- 
iality parameter r = (6 — c)/(a — c) in the range of 0.15 to 
0.25. In the central region, i.e. the inner most 10 per cent, the 
decrease of the axis ratios begins at the same time as at the 
half mass radius, but they are typically larger by 0.1-0.15 (Fig. 
2b). The flattening of the halo region, i.e. the 90 per cent most 
bound particles, is weak and strongly delayed according to the 
enhanced collapse time in the outer region. Additionally, the 
isotropic redistribution of particles gaining energy in the vio- 
lent relaxation phase leads to a more spherical structure of the 
halo. 



Fig. 2b. Temporal evolution of the minor-to— major axis ratio in the 
centre (10% most bounded particles, upper curves after collapse) and 
at the half— mass radius. The softening length was set to e = 0.01 
(dashed, model Al) and e = 0.1 (solid, model Bl) with N = 32768. 
Note here and in the following figures that the timescale of the first 
collapse of the central region is constant and is equal to iff (ro) ta 0.84 
for all models. 
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Fig. 2c. Temporal evolution of the Lagrange radii for a model with 
iV=4096, no softening (+, HARP) and a simulation with e = 0.01 
(solid, GRAPE, model D). 



Fig. 2d. As Fig. 2c, but for A=16384 (GRAPE, model F). 



2.2. A parameter study 

In order to check the influence of the parameters used in N- 
body simulations we selected one reference model (Al) and 
then varied its parameters in a sequence of models. Physically 
we changed the gravitational softening length e (models Bl-3 
and C), the total number N of the particles (models D, F), 
and e and N simultaneously so as to keep the ratio of e to 
the mean interparticle distance constant (model E). In order 
to study the influence of softening on the long-term evolution 
of the quasi-equilibrium configuration formed after ROI, two 
simulations with different softening, but identical configura- 
tions as model F at t = 20 were performed (models G, H). For 
the two models I and J the initial virial coefficient was changed 
from the reference value r\ — 0.04 to r\ — 0.2 in order to com- 
pare our models with the results of Aguilar & Merritt (1990). 
To test the intrinsic errors of our program we varied in mod- 
els A2 and B2 the random number set used to initialize the 
particle configuration, and in model A3 and B3 a significantly 
larger (constant) time step was used (increased by a factor of 
10 as compared to the reference model). For model A3 with the 
small softening length of 0.01 it led to untolerable errors of the 
energy "conservation" (~ 90%), so this model was discarded 
for the data evaluation (see also Sect. 2.2.4). 



For all models with small initial virial coefficient the col- 
lapse of the core proceeded in one free-fall time; Figs. 2c, d 
depict the evolution of Lagrangian mass shells containing the 
indicated fraction of total mass for N = 4096 and N = 16384. 
For comparison data of a run using HARP-2 and a differ- 
ent A'-body integrator (NBODY6 + + , Spurzem 1998) without 
any softening are given. The excellent agreement of both re- 
sults tells us that for the reference model there are no more 
softening-dependent effects in the simulation during the col- 
lapse timescale up to the maximum particle number used in 
our HARP simulations (16384), hence the softened model is 
dynamically equivalent to the "real" discrete system. This how- 
ever, has to be proven again when increasing N further. 

Only after the maximum central density in the collapse 
is reached strong deviations of the shape of the system from 
spherical symmetry (generally triaxiality) occur (as can be seen 
by comparing Figs. 2b and 2c, d). According to Hozumi, Fuji- 
wara & Kan-ya (1996) this could be ascribed to the ROI pro- 
ducing strongly anisotropic velocity dispersions, which then 
give rise to a non-spherical equilibrium system. 

Tables 1-3 give an overview for all models presented here, 
first the model parameters (Table 1), and thereafter some typ- 
ical values describing the flattening (small and intermediate 
axes relative to the major axis, and its time derivative) with 
the timescale t mln defined as the time of minimum c/a (which 
is almost identical to the time of minimum b/a). Note that i m i n 
is larger than the collapse time of the core. 

As a general conclusion from Tables 2 and 3 one could al- 
ready deduce that small softening tends to increase the rate 
of change of c/a after the system has become flattened. At 
a first glance this could be interpreted as a stronger effect of 
two-body relaxation present in these systems. However, we 
will argue in the following that more caution is necessary in 
the interpretation of the data. 

Table 1. Model parameters for the numerical simulations 



model 


N 


s 




At 


comment 


Al 


32768 


0.01 


0.04 


0.0015 


reference model 


A2 


32768 


0.01 


0.04 


0.0015 


different random set 


A3 


32768 


0.01 


0.04 


0.015 


large energy error 


Bl 


32768 


0.1 


0.04 


0.0015 




B2 


32768 


0.1 


0.04 


0.0015 


different random set 


B3 


32768 


0.1 


0.04 


0.015 




C 


32768 


0.02 


0.04 


0.0015 




D 


4096 


0.01 


0.04 


0.0015 




E 


4096 


0.02 


0.04 


0.0015 




F 


16384 


0.01 


0.04 


0.0015 




G 


16384 


0.02 


0.04 


0.0015 




H 


16384 


0.1 


0.04 


0.0015 




I 


32768 


0.01 


0.2 


0.0015 




J 


4096 


0.01 


0.2 


0.0015 





2.2.1. The softening 

For the investigation of the effect of softening we increased e 
by a factor 10 to 0.1 which is a factor of 5 larger than the mean 
central interparticle distance (model Bl). Fig. 2b shows that 
the deviation from sphericity for half of the particles is now 
delayed by a factor 2.5-3 compared to the reference model Al 
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Table 2. Flattening parameters of the 50% most bound parti- 
cles 



model 




c 

a min 


b 

a, min 


c 

a fin 


b 

a fin 


Al 


2.9 


0.36 


0.41 


0.48 


0.57 


A 9 


Q 


U.OO 


u.^ri 


u.oz 


u.oo 


Bl 


6.6 


0.34 


0.39 


0.35 


0.42 


B2 


7.1 


0.37 


0.39 


0.38 


0.42 


B3 


6.6 


0.34 


0.38 


0.38 


0.44 


C 


4.0 


0.34 


0.39 


0.46 


0.57 


D 


3.3 


0.39 


0.42 


0.54 


0.57 


E 


4.3 


0.35 


0.41 


0.50 


0.58 


F 


3.7 


0.35 


0.40 


0.52 


0.57 



The models G and H are omitted because they start with the 
configuration of model F at t = 20. The dynamically hot mod- 
els I, J are omitted since they do not flatten significantly. Model 
A3 is omitted because energy conservation is not sufficient. The 
final time corresponds here to t = 20. 



Fig. 3. Temporal evolution of the minor— to-major axis ratio in the 
centre (10% most bounded particles) and at the half-mass radius. 
The number of particles has been set to N = 32768 (solid, model Al) 
and N = 4096 (dashed, model D) with a softening length e = 0.01. 



Table 3. Flattening parameters of the 10% most bound parti- 
cles 



model 


tmin 




b 

a min 


a fin 


b 

a fin 


Al 


2.9 


0.53 


0.55 


0.62 


0.64 


A2 


4.1 


0.53 


0.56 


0.65 


0.71 


Bl 


6.8 


0.78 


0.80 


0.77 


0.80 


B2 


7.0 


0.76 


0.79 


0.76 


0.79 


B3 


6.5 


0.78 


0.80 


0.77 


0.78 


C 


4.0 


0.54 


0.56 


0.60 


0.63 


D 


2.1 


0.67 


0.72 


0.81 


0.87 


E 


4.3 


0.66 


0.70 


0.67 


0.71 


F 


2.8 


0.60 


0.61 


0.71 


0.77 



Models G,H,I,J,A3 are omitted for the same reasons as in Table 
2. The final time corresponds here to t = 20. 



Fig. 4. Temporal evolution of the minor— to-major axis ratio in the 
centre (10% most bounded particles) and at the half-mass radius. 
The number of particles has been set to N = 32768 (solid, model Al) 
and N = 4096 (dashed, model E) with a softening length e = 0.01 
(model Al) and e = 0.02 (model E), respectively. 



with a smaller softening. The minimum axis ratio c/a reaches 
0.34 &tt = 6.65 for the 50 per cent most bound particles which 
is almost the same minimum value as in model Al. Until the 
end of the simulation c/a increases only very slightly to 0.35. 

Here we find two quite different effects of a variation of 
softening. First, the build-up of non-sphericity needs more time 
if softening is larger. This is valid for both the 10% and the 
50% most bound particles. Second the difference between the 
shape of the 10% and the 50% most bound particles is smaller 
for small softening. Or, in other words, decreasing the softening 
makes the 10% most bound particles less spherical, but the 50% 
most bound particles more spherical. We interpret such results 
as follows: first the initial collapse proceeds to much deeper 
central potential in the case of small softening. So the effect of 
the ROI is stronger, which explains the difference in t m i n of the 
two cases. The term stronger here means that the timescale for 
the build-up of deviations from spherical symmetry is smaller, 
as well as a smaller numerical value of c/a for the 10% most 
bound particles. 



2.2.2. The number of particles 

The effect of a variation of the particle number at constant 
softening is shown in Fig. 3. While until the time i m in differ- 
ences are very small, ROI proceeds to stronger perturbations 
of the spherical shape in the inner shells (10%) for larger N. 

In the models D and E we decreased the total number of 
particles to N = 4096 and performed the calculation with two 
different softening lengths: e = 0.01 (model D) and e = 0.02 
(model E). Until c/a reaches its minimum value at t — 3.0 
the temporal evolution of the axis ratio inside the 50 per cent 
region are very similar to the reference model Al with the 
same softening, but an increased particle number of N = 32768 
(Fig. 3). We conclude that until t m i n the ROI measured at the 
half-mass radius does not depend strongly on the number of 
particles, but only on the depth of central potential reached. 
This supports that the onset of ROI is a collective effect. 

In model E the softening length is a factor 2 larger than in 
model Al and D. This means that the ratio of the softening 
length to the average particle distance is the same for the ref- 
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erence model Al and model E. In the central region the larger 
softening induces a finally increased c/a of 0.7 (Fig. 4). As in 
the previous case larger N generally increases the strength of 
ROI in the inner parts of the system. However, small differences 
in the evolution of the 50% mass shell can be compensated by 
our variation of e ~ l/\fN. This suggests that gravitational 
scattering plays a role here. However, classical two-body re- 
laxation by particle-particle interaction (e.g. Spitzer & Hart, 
1971) is too slow to explain the observed compensation. Thus, 
the scattering of particles with the massive and much more 
concentrated core of the system as a whole might cause the 
relaxation here. This point will be discussed in Sect. 3 in more 
detail. We are grateful to T. Tsuchiya (pers. communication) 
to point this out to us as a possible mechanism. 

2.2.3. The virial coefficient 



For several models we performed the simulations with other 
sets of initial configurations (e.g. models A2 and B2). Gener- 
ally, the deviations of the axis ratios are less than 5 per cent, 
until dynamical equilibrium has been established, i.e. until 
iw5. However, the relaxation timescale seems to vary slightly 
which gives an intrinsic error of c/a of 0.04 at the end of most 
of the simulations (t = 20). Thus, our results do not critically 
depend either on the chosen timestep or on the randomly cho- 
sen initial particle configuration. 

2.3. The long-term experiments 



Fig. 5. Temporal evolution of the minor-to— major (c/a) and inter- 
mediate-to— major (b/a) axis ratio for the 50% most bounded parti- 
cles. The initial virial coefficient has been set to n v ir = 0.04 (solid, 
model Al) and r) v i T = 0.2 (dashed, model I). Both models were 
calculated with N = 32768 and e = 0.01. 

We studied the influence of the initial distribution in veloc- 
ity space by starting with a dynamically "warm" configura- 
tion for the spatial configurations of the models Al and D, i.e. 
a virial coefficient rj v ir = 0.2 (models I (N = 32768) and J 
(N = 4096)). In agreement with Aguilar & Merritt (1990) and 
Boily et al. (1998) the memory to the initial conditions is pre- 
served and the system remains almost spherical throughout the 
simulations when starting from spherical initial configurations 
(Fig. 5) . Within statistical errors the models do not depend on 
the total number of particles. 

2.2.4. Miscellaneous 

In a set of models we investigated the influence of numerical 
effects introduced by the chosen timestep (models A3 and B3) 
and the random realization of the initial particle distribution 
(models A2 and B2). If the timestep is increased by an order 
of magnitude (model A3) the energy is only well conserved 
(better than 1.5 per cent) for a model with a large softening 
of £ = 0.1 (H), whereas the small softening of the reference 
model introduces a very large integration error at the moment 
of strongest compression. If we perform the simulation of model 
Al with a doubled timestep, the results are almost identical to 
the reference model. 



Fig. 6. Long-term evolution of the intermediate-to-major axis ratio 
(solid) and the minor-to-major axis ratio (dashed) at the half-mass 
radius for a simulation with N = 16384 particles. The softening 
length was set to e = 0.01 (model F). 



Fig. 7. Long-term evolution of the minor-to— major axis ratio at 
the half-mass radius for simulations with N = 4096 (solid, model 
D) and N = 16384 (dashed, model F) particles. The softening length 
was set to e = 0.01. The straight lines correspond to a sphcrization 
of the system in one relaxation time (the latter is already corrected 
for softening). 



In some cases we followed the evolution of the system for a 
much longer time, until it practically relaxed again to a spher- 
ically symmetric configuration (models Al, C, D, E, F, G and 
H). The evolution of shape shows three different phases (e.g. 
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Fig. 8. Long— term evolution of the (logarithmic) deviation from 
sphericity at the half— mass radius for simulations with N = 4096 
(solid, model D) and N = 16384 (dashed, model F) particles. The 
softening length was set to e = 0.01. 



Fig. 9. Long— term evolution of the (logarithmic) deviation from 
sphericity at the half-mass radius for simulations with a softening 
length e = 0.01 (solid, model F) and e = 0.02 (dashed, model G) 
particles. The simulation was performed with N = 16384 particles. 



Figs. 6 and 11): first after the collapse the ROI leads to a tri- 
axial configuration. In the second stage a fast readjustment 
of shape occurs during about 20 crossing times. Thereafter a 
steady relaxation follows until the system becomes axisymmet- 
ric (b/a ^ 1, but still c/a 1) and finally spherical (Fig. 6). 
Lateron the system does not change its shape (Fig. 7) . 

The timescale for spherization is given remarkably exact by 
the two-body relaxation time t T i x of the N-body system. This 
is demonstrated in Fig. 7 by the temporal evolution of the 
minor-to-major axis ratio for two simulations with N = 4096 
and N = 16384 and a softening length of e = 0.01. The straight 
lines are no fits, but the slopes one gets by the assumption 
that spherization acts on a relaxation timescale. It should be 
noted that the good agreement can only be achieved when the 
Keplerian relaxation time is corrected for softening (for details 
see Theis (1998)). This gives a correction, i.e. an increase of the 
relaxation times, by approximately a factor of 3 for N = 16384. 

The deviation from sphericity, S(t) = 1 — c/a(t), can be 
fitted by S(t) « 5(0) exp(— £/t r ix) as the extended linear regime 
in Fig. 8 shows. Both, the dependence of the axis ratio evolu- 
tion on the total number of particles (Fig. 8) and the softening 
length (Fig. 9 and 10) are in agreement with a spherization 
within a relaxation timescale of the softened two-body poten- 
tial. The slow relaxation in model H (e = 0.1) is due to the 
small ratio of maximum impact parameter to softening length 
(e is twice the core radius!) which gives a strongly enhanced 
relaxation time. 

Shortly after the initial collapse at i m in, however, the "re- 
laxation" of shape of our models is much faster than one could 
expect from the two-body relaxation argument. This can be 
seen clearly in Fig. 11, which shows the long-term evolution 
of our reference model Al (N = 32768, £ = 0.01). We denote 
this phase of fast shape change by the term "violent shape re- 
laxation". It is charcterized by a substantially larger growth 
rate of the axis ratio as compared with expectations from two- 
body relaxation. Violent shape relaxation sets in after the ROI 
has established a triaxial configuration. Core radius and core 
mass do not vary strongly during violent shape relaxation (see 
Fig. 2a). Independent of taking 10% or 50% of the particles, 
the axis ratio growth rate during violent shape relaxation is 
constant for about 10-20 crossing times; thereafter a smaller 
growth rate obtained by two-body relaxation is following. Note 
that all variations of shape seen in Figs. 2-5 and discussed in 
subsection 2.2 are still in the violent phase and do not rep- 
resent the long-term shape changes from two-body relaxation, 
which can be seen in the Figs. 6-11 for the long-term evolution. 



3. Discussion and Summary 



Fig. 10. Long— term evolution of the (logarithmic) deviation from 
sphericity at the half-mass radius for simulations with a softening 
length e = 0.1 (solid, model H) and e = 0.02 (dashed, model G) 
particles. The simulations were performed with TV = 16384 particles. 



We performed a set of A-body simulations for the collapse 
of a dynamically cold Plummer sphere and studied the oc- 
currence of radial orbit instability (ROI) and the evolution of 
shape. Especially, we focussed on the influence of the standard 
N-body parameters like the total number of particles or the 
softening length. Our simulations showed that the growth rate 
of the ROI is strongly changed by the softening at constant 
particle number and almost unaffected by a variation of N at 
constant softening. For our standard models with e = 0.01 and 
N — 16384 (N — 4096), however, we could show by comparison 
with high-accuracy zero softening models, that no more arti- 
ficial softening-dependent effects are present. But our results 
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Fig. 11. Long-term evolution of the minor-to— major axis ratio in 
the centre (10% most bounded particles, dashed) and at the half- 
mass radius (solid) for model Al (N=32768, e = 0.01). 



Fig. 12. Temporal evolution of the distance of the density centre 
from the centre of mass (solid) in correlation with Lagrangian radii 
(0.1%, 20%, 50%, dashed) for model Al (N=32768, e = 0.01). 



Fig. 13. Long-term evolution of the minor— to-major axis ratio at 
the half-mass radius for simulations with a softening length e = 0.01 
(solid, model Al) and e = 0.02 (dashed, model C). The simulations 
were performed with N = 32768 particles. 



well as the undesired two-body relaxation. Such findings are 
consistent with recent arguments given by Steinmetz & White 
(1997), Moore et al. (1997), and Fukushige & Makino (1997), in 
the context of iV-body simulations for cosmological structure 
formation. 

Generally, according to our simulations, the cold collapse 
with ROI and the subsequent evolution could be subdivided 
into three phases. First there is an initial strong collapse and 
build-up of triaxiality via ROI, characterized by a time i m i n at 
which a minimum value of the principal axis ratios c/a and b/a 
is obtained. Here our results agree with the findings of Hozumi, 
Fujiwara & Kan-ya (1996) insofar, that the system remains 
spherical until maximum collapse t co ii < t m in; it is only at that 
stage that strong anisotropies of the velocity dispersions create 
a triaxial configuration. This is also consistent with arguments 
by Bettwieser & Spurzem (1986) and Hensler et al. (1995), 
who state that homologously collapsing systems have to stay 
isotropic. 

Second, a phase of rapid, violent shape relaxation occurs. 
It is much faster than by any estimate due to two-body re- 
laxation. The mechanism of such violent shape relaxation is 
still unclear. We observe, however, strong fluctuations of the 
density centre of the system just during that phase, which are 
strongly correlated with oscillations of the Lagrangian radii of 
up to 20% of the total mass (see Fig. 12, here Lagrangian radii 
have not been measured with respect to the density centre but 
with respect to the centre of mass of the system). Even the 
half-mass radius varies due to these quasi-oscillations of the 
core whereas the 0.1% Lagrange radius is completely inside the 
core and, thus, not affected by the core's motion. We conclude 
that the rapid initial collapse induces strong motions of the 
density centre, which decay on a timescale of some ten cross- 
ing times. However, we saw that the core is already formed 
after two crossing times and does not change any more. This 
means, that a core fixed in size moves in the central region, by 
this inducing strong potential fluctuations which give rise to a 
relaxation which is added to the unavoidable two-body relax- 
ation. When the core settles down these large-scale potential 
fluctuations disappear and the additional relaxation does not 
operate any more by this defining the last stage of shape evo- 
lution. A comparison of two long-term simulations with differ- 
ent softening lengths strengthens that point, insofar as it seems 
that two-body relaxation does not change the violent shape re- 
laxation (Fig. 13): Whereas a difference in the relaxation rate 
in the last stage after t = 20 is obvious, the shape evolution is 
qualitatively almost unaffected during the violent shape relax- 
ation stage. However, the relatively short timescale of violent 
shape relaxation and the scatter in the measured c/a-values 
makes a quantitative comparison impossible. It is interesting 
to compare with the similar problem of core oscillations in 
quasi-statically evolving spherical star cluster undergoing core 
collapse (Makino & Sugimoto 1987, Heggie, Inagaki & McMil- 
lan 1994, Spurzem & Aarseth 1996), which seem to occur on a 
similar timescale. 

Finally, a third phase is observed, which is characterized by 
a secular drift of the system via axisymmetry towards spherical 
symmetry. It occurs on a time scale consistent with that of 
two-body relaxation, provided the proper softening is taken 
into account (Theis 1998). 



show that for a given particle number inappropriate softening 
could yield wrong results by changing the collective effects as 
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